####################################################################
## author:    Robert A. Huber
## contact:   robert.huber@ir.gess.ethz.ch
## file name: pc_uk_sem.R
## Context:   Populism and Climate Sceptism, individuals from BES
## started:   2019-05-19
## Summary:   runs SEM
######################################################################

library(lavaan)

#Prepare variables ####

df_bes$pop1sem <- ifelse(df_bes$populism1 == "Don't know", NA,
                         as.numeric(df_bes$populism1))

df_bes$pop2sem <- ifelse(df_bes$populism2 == "Don't know", NA,
                         as.numeric(df_bes$populism2))

df_bes$pop4sem <- ifelse(df_bes$populism4 == "Don't know", NA,
                         as.numeric(df_bes$populism4))

df_bes$pop5sem <- ifelse(df_bes$populism5 == "Don't know", NA,
                         as.numeric(df_bes$populism5))

df_bes$pop6sem <- ifelse(df_bes$populism6 == "Don't know", NA,
                         as.numeric(df_bes$populism6))

df_bes$auth1sem <- ifelse(df_bes$auth1 == "Don't know", NA,
                               as.numeric(df_bes$auth1)-1)

df_bes$auth2sem <- ifelse(df_bes$auth2 == "Don't know", NA,
                               abs(as.numeric(df_bes$auth2)-2))

df_bes$auth3sem <- ifelse(df_bes$auth3 == "Don't know", NA,
                               abs(as.numeric(df_bes$auth3)-1))

df_bes$auth4sem <- ifelse(df_bes$auth4 == "Don't know", NA,
                               abs(as.numeric(df_bes$auth4)-1))

df_bes$efficacy1sem <- ifelse(df_bes$efficacyPolCare == "Don't know", NA,
                                             as.numeric(df_bes$efficacyPolCare))

df_bes$efficacy2sem <- ifelse(df_bes$efficacyNoMatter == "Don't know", NA,
                                              as.numeric(df_bes$efficacyNoMatter))

df_bes$efficacy3sem <- ifelse(df_bes$efficacyNotUnderstand == "Don't know", NA,
                                                   as.numeric(df_bes$efficacyNotUnderstand))

df_bes$efficacy4sem <- ifelse(df_bes$efficacyUnderstand == "Don't know", NA,
                                                -as.numeric(df_bes$efficacyUnderstand))

df_bes$efficacy5sem <- ifelse(df_bes$efficacyTooMuchEffort == "Don't know", NA,
                                                   as.numeric(df_bes$efficacyTooMuchEffort))

#+ auth + gender + inc + age  + edu_high + riskTaking + polAttention + satDem
#+ efficacy + genTrust + econPersonalRetro +  econGenRetro
#+ partyId_ukip + partyId_con + partyId_lab + partyId_lib + partyId_gre + partyId_oth

df_bes[,c("climateChange_ord", "enviroProtection_ord")] <- lapply(df_bes[,c("climateChange_ord", "enviroProtection_ord")], ordered)


### BAseline

modelCFA <- '
   # latent variables
     pop =~ pop1sem + pop2sem + pop4sem + pop5sem + pop6sem
     auth =~  auth1sem + auth2sem + auth3sem + auth4sem
     eff =~ efficacy1sem + efficacy2sem + efficacy3sem + efficacy4sem + efficacy5sem
   # residual covariances
      efficacy3sem ~~ efficacy4sem
      pop2sem ~~      pop4sem
      efficacy3sem ~~ efficacy5sem
      efficacy4sem ~~ efficacy5sem
      pop1sem ~~      pop2sem
'

fit <- cfa(modelCFA, data= df_bes, meanstructure = T)
summary(fit, standardized = F, fit.measures = T, rsquare = T)
modindices(fit, sort = T, minimum.value = 3.85)

fitmeasures(fit)[c("chisq", "df", "cfi", "tli", "rmsea", "srmr")]

## Climate Scepticism ####

df_bes$climateChange_ord <- factor(df_bes$climateChange_ord,
                                      ordered = T)

modelSEMCC <- '
   # latent variables
     pop =~ pop1sem + pop2sem + pop4sem + pop5sem + pop6sem
     auth =~  auth1sem + auth2sem + auth3sem + auth4sem
     eff =~ efficacy1sem + efficacy2sem + efficacy3sem + efficacy4sem + efficacy5sem
   # regressions
     climateChange_ord ~ lr + pop + auth + eff + age + gender + inc + econPersonalRetro + econGenRetro + satDem
   # residual covariances
      efficacy3sem ~~ efficacy4sem
      pop2sem ~~      pop4sem
      efficacy3sem ~~ efficacy5sem
      efficacy4sem ~~ efficacy5sem
      pop1sem ~~      pop2sem
'

fit <- sem(modelSEMCC, data = df_bes, sampling.weights = df_bes$wt_full_w7)
summary(fit, standardized=TRUE, fit.measures = T)
mi <- modindices(fit, minimum.value = 3.85, sort=T)
parameterEstimates(fit)[15:20,]
fitmeasures(fit)[c("chisq", "df", "cfi", "tli", "rmsea", "srmr")]

## Environmental Protection ####

df_bes$enviroProtection_ord <- factor(df_bes$enviroProtection_ord,
                                      levels = c("Gone too far", "About right", "Not enough"), ordered = T)

modelSEMEP <- '
   # latent variables
     pop =~ pop1sem + pop2sem + pop4sem + pop5sem + pop6sem
     auth =~  auth1sem + auth2sem + auth3sem + auth4sem
     eff =~ efficacy1sem + efficacy2sem + efficacy3sem + efficacy4sem + efficacy5sem
   # regressions
     enviroProtection_ord ~ lr + pop + auth + eff + age + gender
   # residual covariances
           efficacy3sem ~~ efficacy4sem
      pop2sem ~~      pop4sem
      efficacy3sem ~~ efficacy5sem
      efficacy4sem ~~ efficacy5sem
      pop1sem ~~      pop2sem
'

fit <- sem(modelSEMEP, data = df_bes, sampling.weights = df_bes$wt_full_w7)
summary(fit, standardized=TRUE, fit.measures = T)
mi <- modindices(fit, minimum.value = 3.85, sort=T)
parameterEstimates(fit)[15:20,]
fitmeasures(fit)[c("chisq", "df", "cfi", "tli", "rmsea", "srmr")]

## Environment vs Growth ####

modelSEMEG <- '
   # latent variables
     pop =~ pop1sem + pop2sem + pop4sem + pop5sem + pop6sem
     auth =~  auth1sem + auth2sem + auth3sem + auth4sem
     eff =~ efficacy1sem + efficacy2sem + efficacy3sem + efficacy4sem + efficacy5sem
   # regressions
     enviroGrowth ~ lr + pop + auth + eff + age + gender
   # residual covariances
      efficacy3sem ~~ efficacy4sem
      pop2sem ~~      pop4sem
      efficacy3sem ~~ efficacy5sem
      efficacy4sem ~~ efficacy5sem
      pop1sem ~~      pop2sem
'

fit <- sem(modelSEMEG, data = df_bes, sampling.weights = df_bes$wt_full_w7)
summary(fit, standardized=TRUE, fit.measures = T)
mi <- modindices(fit, minimum.value = 3.85, sort=T)
parameterEstimates(fit)[15:20,]
fitmeasures(fit)[c("chisq", "df", "cfi", "tli", "rmsea", "srmr")]